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ABSTRACT 

The measured properties of stellar oscillations can provide powerful constraints on 
the internal structure and composition of stars. To begin this process, oscillation fre- 
quencies must be extracted from the observational data, typically time series of the 
star's brightness or radial velocity. In this paper, a probabilistic model is introduced 
for inferring the frequencies and amplitudes of stellar oscillation modes from data, 
assuming that there is some periodic character to the oscillations, but that they may 
not be exactly sinusoidal. Effectively we fit damped oscillations to the time series, 
and hence the mode lifetime is also recovered. While this approach is computationally 
demanding for large time series (> 1500 points), it should at least allow improved 
analysis of observations of solar-like oscillations in subgiant and red giant stars, as 
well as sparse observations of semiregular stars, where the number of points in the 
time series is often low. The method is demonstrated on simulated data and then 
applied to radial velocity measurements of the red giant star ^ Hydrae, yielding a 
mode lifetime between 0.41 and 2.65 days with 95% posterior probability. The large 
frequency separation between modes is ambiguous, however we argue that the most 
plausible value is 6.3 fiHz, based on the radial velocity data and the star's position in 
the HR diagram. 

Key words: stars: oscillations — methods: statistical — stars: individual: £, Hydrae 



1 INTRODUCTION 

The study of stellar oscillations provides a powerful probe 
of the physical properties of stars. In particular, knowledge 
of the frequencies of many eigenmodes of a star can sig- 
nific antly constrain its internal structure and composition 
(e.g. iHoudekllioOTi ). In practice, these frequencies are in- 
ferred from time series data of the star's radial velocity or 
intensity, which is analysed in order to determine the fre- 
quenci es of the oscillation mode s that contributed to the 
signal (|Bedding fc Kield sen 2007). The amount and quality 
of data has increased spectacularly over the past few years, 
mostly due to advances in instrumentation t hat were primar- 
ily in tended for extrasolar planet searches (|Marcv fc Butlen 
1 19921 ). Despite these advances, oscillation data on stars other 
than the Sun is still much more sparse and noisy, for ob- 
vious reasons. Most of this data is analysed with Fourier 
p ow er spectrum methods inherited from helioseismolog y 
(e.g. iToutain fc Froehlich|[l99i : Ijimenez-Reves et al.|[20oi ). 
Recently, there has been growing interest in examining the 
fundamentals of data analysis techniques, and attempts to 
improve on these classical techniques have yielded modest, 
but non-negligible improvements to our ability to make use 
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of time series data (e.g. iFletcher et al.l[2006l: iBrewer et al.l 



2007'; iRe gulo fc Roca Corted '200?': ' Stahn fc GizonI \20m . 
Grubcrb auer. Kallinger. fc W eiss 2008]). Separately, there 
has been a steady grow th in interest in Bayesian Inference 
jSivia and SkiUineH2006i ) as the most consistent and natural 
way to model uncertainties. Thus, the approach described 
in this paper is Bayesian. 

The idea behind Bayesian methods is to describe our 
knowledge by a probability distribution over the space of 
possible solutions we are considering. This probability dis- 
tribution then gets updated to take into account the infor- 
mation contained in the data, in this case the time series 
data In asteroseismology, the possible solutions we 

consider are all possible values for the parameters of in- 
terest: the frequencies {i/} of the oscillation modes, their 
amplitudes {A} and phases {4>}, and the total number of 
modes, m. For brevity, we drop the braces hereafter; A, 
V and (\> now stand for arrays of amplitudes, frequencies 
and phases respectively. Any additional parameters (mode 
lifetime, for example) are denoted collectively by d. Before 
taking into account the data, we assign a prior distribu- 
tion p(6,m,v, A, (j)). We also probabilistically model the pre- 
dictions for what data y we expect to observe as a func- 
tion of the parameters: p{ii\6,m,i',A,ij)), called the sam- 
pling distribution. Sometimes it is possible to marginalise 
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out the amplitude and phase parameters, by integrating over 
them. This helps the computational search, because the al- 
gorithms only need to find good values for the frequencies, 
not the a mplitudes and phase s. This was done in a previ- 
ous paper l|Brewer et ahlfioOTl . hereafter B07) by replacing 
the amplitudes and phases with sine-amplitudes and cosine- 
amplitudes; however, it is not possible for the model we in- 
troduce in this paper. Throughout this paper, we will be 
fitting both the frequencies and the amplitudes. It turns out 
(Section 12. ip that we do not need phases in our model, so 
we will drop the phases hereafter. 

Once the prior distribution and the sampling distribu- 
tion have been specified, we have defined a joint probability 
distribution for the parameters and the data: 

p{e, m, V, A, y) = p{e, m, u, A)p{y\6, m, u, A) (1) 

This probability distribution describes what we know about 
the parameters and the data before we observe the actual 
data. Once we learn the actual values of the data Vq^^^, we 
update our probability distribution by deleting all hypothe- 
ses in the (6', m, v, A, y) space that are now known to be false; 
i.e. we restrict our attention to the slice y — j/gjjg. This gives 
the posterior distribution for the parameters given the data, 
describing our knowledge after updating to include the effect 
of the data. This is expressed by Bayes's theorem: 

p{e,m,u,A\y = yobs) p{S,m, u, A)p{y\e,m, u, A)\y^^^ (2) 

The second factor in Equation [5] is the sampling distribu- 
tion (probability distribution for the data as a function of 
the parameters) with the data fixed at the observed values; 
thus it is a function of the parameters only. Once the data 
have been fixed, it is commonly referred to as the likelihood 
function. 

Many existing methods, including the Bayesian Analy- 
sis of B07, rely on the assumption that the observed time se- 
ries is composed of a sum of sinusoidal signals, plus Gaussian 
noise - basically, this was our choice for the sampling distri- 
bution p{y\9,m, ly, A). Interestingly, the periodograirQ can 
be proven to be a sufficient statistic from this same assump- 
tion, plus the assumption that the time series has complete 
phase coverage. What this means is that if we intend to infer 
the frequencies of the modes, no informat ion is lost by re- 
ducing the tim e series to the periodogram (|Bretthorstll 19881 : 
lGregorvl[200ll ). as long as the signal is purely sinusoidal and 
the time series has no significant gaps. Successful methods 
have been developed to infer frequencie s and mode lifetimes 
by fi tting to the power spectrum (e.g. lAppourchaux et al.l 
l2008l ). However, the presence of stochastically excited modes 
and gaps in the data both break the assumptions required 
for the periodogram to be a sufficient statistic. Hence, ide- 
ally, when either of these conditions do not hold, we should 
work with the raw time series data. 

Of course, it may be the case that using the power spec- 
trum discards an insignificant amount of information, and 
the gain in convenience of the power spectrum far outweights 
such theoretical concerns. This is certainly the case with so- 
lar data, and stellar data with good coverage and a long 

^ Throughout this paper, the terms "periodogram" and "power 
spectrum" will be used interchangeably. It is also common to plot 
the square root of the periodogram, which is sometimes called the 
"amplitude spectrum" . 



mode lifetime. It may be more generally true; however, this 
question is beyond the scope of this paper. In this paper 
we describe a new method to analyse time series data, tak- 
ing into account the fact that the predicted signal from an 
oscillation mode is quasi-sinusoidal, and that the data may 
contain gaps, removing any concerns about information loss 
due to pre-processing via taking the power spectrum. This 
is done by making the choice for p{y\9,m, v,A) as realistic 
as possible. 



2 SOLAR-LIKE OSCILLATIONS 

In the following, we will consider a single stochastically ex- 
cited mode, and obtain a model for the sampling distribution 
p{y\6,m,v, A) when m = 1. Subsequently, we will expand 
this to include multiple oscillation modes and observational 
errors, in Section [2.21 Solar-like oscillations are oscillations 
in main sequence or red giant stars that are continually 
damped and re-excited by turbulent convection, and there- 
fore do not produce purely sinusoidal signals. For solar-like 
oscillations, the signal due to a single mode is modelled as 
a damped and stochastically excited oscillator with a driving 
force fit): 

where t is a damping timescale (the factor of two is intro- 
duced such that solutions to Equation [3] without the driving 
force decay with an e-folding time of r), and 13 is an ampli- 
tude constant for the driving force f{i). If f{i) is specified 
and the initial conditions j/(0) and y{Q) are known, then 
Equation [3] has a unique solution. However, as the term 
stochastic excitation suggests, f{t) is only specified proba- 
bilistically. Throughout this paper, we will assume that fit) 
is unit variance white noise, so f{t) at any time t comes 
from a standard Gaussian distribution with mean and 
standard deviation 1. /(fi) and /(t2) are independent for 
all distinct times ti 7^ i2. The white noise probability distri- 
bution that we assigned to fit) is an example of a Gaus- 
sian Process distribution. In general, a Gaussian Process 
is a p robability distribution over a space of poss ible func- 
tions (jRasmussen fc WiUiams|[200^ : iMacKavlliooal ) ; however 
a more general Gaussian Process may differ from white noise 
because the function value at different times may be corre- 
lated. 

Since Equation [3] is a linear ordinary differential equa- 
tion, and fit) is assigned a Gaussian process distribution, 
we have implicitly also assigned a Gaussian process distribu- 
tion for the value of the oscillating signal y{t) at all times t. 
Whereas /(ti) and /(t2) are independent (for ti 7^ t2), the 
value of the oscillation signal at any two times, j/(ti) and 
yit2), are correlated, with the covariance function defined 
as 

C{t,,t,) = {yit,)yit,)) - {y{t,)) {yitj)) 

= (y(iOs/(ii)> (4) 

where the expectatation value (mean) of y, {y{t)), has been 
set to zero for all time. A typical signal obtained by solving 
Equation 13] is shown in Figure [T] Clearly, the signal value at 
any given time is strongly correlated with the signal value 
at a time one period later, and anticorrelated with the value 
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half a period later. However, the correlation is weaker than 
would be the case if the signal was purely sinusoidal. The 
entire signal displayed in Figure[T]can be regarded as a single 
point sampled from a Gaussian Process distribution with 
mean function zero and a particular covariance function. 
Equation 13 is capable of producing solutions with any initial 
phase, which is why we do not need phases in our model. 

A useful property of Gaussian Processes is that the 
joint distribution for the function evaluated at a finite set 
of points (i.e. the signal y{t) evaluated at the observation 
times) is a multivariate Gaussian, with covariance matrix 
given by the covariance function evaluated at the relevant 
points. Hence, the joint distribution for the value of the os- 
cillation signal at a discrete set of A*' times {ti, t2, tjv} 
is: 

p{y\A'^,T)= \ ^=exp (-^y'^C'V) (5) 

V(27r)^detC ^2 J 

where y = j/(f2), 2/(tjv)}, the expectation values 

(means) of all of the y's are zero, and C is a covariance ma- 
trix that implicitly depends on A, v and r. Equation[5]is the 
first step in the construction of a realistic p(i/|6, m, v, A): it 
would suit perfectly if we had noise-free data containing one 
mode, and if we knew how C depended on A, v and r. The 
dependence of C on y4, i/ and r is addressed in Section [2. II 
while the generalisation to noisy data with multiple modes 
occurs in Section [2.21 

Throughout this paper, no results depend on the choice 
of the origin for only relative times matter. In this case, 
the covariance function is said to be stationary. This im- 
plies that the covariance of y{ti) and y{t2) depends only on 
the difference between ti and t2, and not on their absolute 
values: 

Citi,t2) = C(ti,t2) = Ci\t2-ti\) = C{At) (6) 

where the symbol C has been used to denote the covariance 
function, whether it takes one argument or two. This result 
is used to evaluate the elements of the covariance matrix C. 

2.1 Use of Simulations to Determine the 
Covariance Function 

In this section, we aim to find exactly how the covariance 
function for the signal depends on A, v and r, for a single 
mode. The dependence on A is trivial: if y{t) comes from 
a Gaussian Process distribution with mean function zero 
and covariance function C{At), then A x y{t) comes from 
a Gaussian Process distribution with mean zero and covari- 
ance function A'^C{At). Note that these quasi-sinusoidal sig- 
nals do not have a strict amplitude like sinusoidal signals do, 
the amplitude A is really just the expected standard devia- 
tion of the oscillation signal. 

To investigate how the covariance function of a stochas- 
tically excited oscillation signal depends on frequency and 
mode lifetime. Equation [3] was solved numerically using a 
fourth order Runge-Kutta algorithm with a timestep much 
smaller than the natural period of the oscillations. From a 
very long simulation, the covariance function for solutions of 
Equation [3] was estimated by taking random pairs of times 
ti and t2, and plotting the average value of y{t{)y(t2) as 
a function of /Si = \t2 ~ t-i\. A short section of the simu- 
lated time series is plotted in Figure [1] which clearly shows 



that while there is some periodic nature to the oscillations, 
the varying amplitude and phase changes would frustrate 
any simple modelling of the signal as sinusoidal waves - too 
many frequencies will be required in order to fit the data 
(B07). 

The estimated covariance function of the signal due to 
a single mode is shown in Figure [2] It can be accurately 
modelled by a cosine curve multiplied by an exponential 
decay: 

C(At) = A^x exp (-|At|/r') cos (27r!/Af) (7) 

where the decay timescale in the covariance function is em- 
pirically found to agree with r to within 5 per cent (Fig- 
ure 0), and in practice we will take it as being equal. In 
the absence of stochastic excitation and damping, the co- 
variance function would be just a cosine function, with fre- 
quency equal to the oscillation frequency. Hence, the effect 
of stochastic excitations and damping can be parameterised 
by the single parameter r and its effect is to put an exponen- 
tial decay factor into the covariance function for the signal. 
In the next subsection, this result is generalised to include 
multiple modes and observational errors. 

2.2 Addition of Independent Gaussian Processes 

Suppose there are two functions of time (for instance, the 
signal from two modes), x{t) and y{t) and our knowledge 
of these functions is described by independent stationary 
Gaussian Processes for each: with mean zero and covariance 
functions Ci(At) and Cy{At) respectively. If we are inter- 
ested in the sum 

zit) = xit)+y{t) (8) 

then the sum is also a Gaussian process with mean zero and 
covariance 

C,(At) = {z{t)z{t + At)) 
= {{x{t) + y{t)) (x(t + At) + y{t + At))) 
= {x{t)x{t + At)) + {y{t)y{t + At)) 
+ {y{t)x{t + At)) + {x{t)yit + At)) (9) 

Since x and y are independent, the expectations of the last 
two terms are zero. Hence 

C,{At) = {x{t)x{t + At)) + {y{t)y{t + At)) 

^ C4At) + Cy{At) (10) 

Therefore, the covariance function for the sum of two inde- 
pendent Gaussian Processes is the sum of their individual co- 
variance functions. This result can easily be extended to any 
number of Gaussian processes. Thus, if we are testing the 
hypothesis that there are many modes with various frequen- 
cies and amplitudes, and that there is also Gaussian noise in 
the data, the relevant covariance matrix (Equations|4]and[5ll 
is the sum of the covariance matrix for each mode (obtained 
from Equation [7)l and a diagonal covariance matrix for the 
noise, with the given measurement uncertainties used as the 
values for the noise standard deviation. In addition to these 
components, we included an "extra noise" signal to account 
for unmodelled errors, misquoted error bars, or correlated 
noise due to stellar effects that are not the oscillations of in- 
terest. The extra noise signal has an unknown standard de- 
viation parameter o-g^tra ^^'^ ^ correlation timescale To- for 
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Figure 1. Simulated signal from a single damped, stochastically excited oscillation mode: a solution to Equation |3] The frequency is 1 
time unit, and the damping timescale (mode lifetime) r is 10 time units. 
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Figure 2. The correlation function for solutions to the damped, stochastically driven oscillator problem, estimated from numerical 
simulations. The correlation is the same as the covariance function but with the variance scaled out. The top panel shows the result 
from simulations with an input frequency of 1 time unit and damping timescale (mode lifetime) of 10 time units. When parameterising 
the covariance function by C(At) = X exp (— | A4|/r') cos {2-KvHt), the best fit value of r' is empirically found to be very close to the 
value of T that was used to produce the data, so r' « r. The bottom panel shows the residuals after subtracting the fitted exponentially 
decaying cosine curve. 



its exponentially decaying covariance function. Thus, o-g-^tra 
and To- are additional parameters to be estimated from the 
data. 



3 BAYESIAN INFERENCE 

We measure the signal at a discrete set of times 
{fi, in} with additive Gaussian noise of standard de- 
viation IV'^^ + '^^xtra'V'^a+'^extra' 



V'^" + '^extra}' 



where 0-1,0-2... are the reported error bars on the obser- 
vations. The probability distribution for the total data set 
given all of the parameters (number of modes, their fre- 
quencies and amplitudes, the extra noise and its timescale) 
is Gaussian with covariance matrix formed by the sum of 
the covariance functions for each component (Section \2.2\ . 
Thus, we have now constructed the sampling distribution 
- the probability distribution for the observed data given 
the parameters of interest. The sampling distribution is the 
same as Equation (5] but where the covariance matrix C is 
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the sum of the covariance matrices for each mode, the diag- 
onal covariance matrix of the measurement errors, and the 
covariance matrix for the extra noise term. The posterior 
probability distribution for the parameters of interest given 
the data is then given by Bayes's theorem (Equation [2]), 
where 9 — {r, cTgxtrai """o"}- i-^- the mode hfetime, extra noise 
standard deviation, and the correlation timescale for the ex- 
tra noise term. 

The previous sections described the sampling distribu- 
tion, and thus the likelihood function, the second term in 
Equation (2] in terms of Gaussian Processes. Now we must 
assign prior distributions for all of the parameters, i.e. the 
first term in Equation [2] For simplicity, we chose the priors 
for all of the parameters to be independent of each other. In 
principle, this could be improved; for example, the expected 
amplitude of a mode is not the same at all frequencies. 

The prior for the number of modes, m, was a uniform 
probability distribution ranging from 1 to a user-specified 
maximum number, which we took to be 200. The prior for 
the frequencies u was a uniform distribution between a user- 
specified lower and upper limit: for the data sets discussed 
in this paper, these limits were and 200 /iHz. The prior for 
the amplitudes A was chosen to be an exponential distribu- 
tion with unknown mean /j,, which effectively becomes yet 
another parameter to be inferred from the data. The priors 
for fi, and the remaining parameters r, cTgxtra ^^'^ '''o-' 
positive parameters, were chosen to be scale-invariant pri- 
ors of the form p{x) oc 1/x between generous upper and 
lower limits. These priors correspond to uniform priors for 
the logarithm of the quantities, and is appropriate for pos- 
itive parameters with unknown order of magnitude. Since 
the specification of the priors introduced an extra parame- 
ter fi to be inferred, the additional parameter vector is now 
extended to include fi. Thus, 9 = {r, o-g^trai '''o-, A*}- 



4 MARKOV CHAIN MONTE CARLO 

The posterior distribution can be effectively sampled using 
Markov Chain Monte Carlo (MCMC). I n our impl ement a- 
tion, we used the Metropolis algorithm (|Neal|[l993l l. Start- 
ing from a model with a single mode of arbitrary frequency 
and amplitude, and typical values for the additional param- 
eters 9 — {r, o-gxtrai '''o'l m}) propose to either add a mode 
(with its frequency and amplitude chosen from the prior), 
remove a mode, move a mode's frequency or amplitude, or 
shift the value of one of the additional parameters such as 
mode lifetime r. Then, the proposed change is accepted with 
a probability that depends on the relative likelihoods and 
prior probabilities of the current and the proposed model. 
Steps to models with higher posterior probability are always 
accepted, steps to models with a lower posterior probabil- 
ity are accepted with a probability given by the ratio of the 
posterior probability of the proposed model to the proba- 
bility of the current one. If a proposed change to the model 
is rejected, the next model in the sequence is the same as 
the previous one. When this algorithm runs, the output of 
the code is a random sequence of models (sets of frequen- 
cies and amplitudes), each possibly slightly different from 
the last, where the diversity amongst the models is indica- 
tive of the uncertainty of any inference. To save memory, a 
subset of effectively independent models from this sequence 



may be used for any subsequent calculations. This is called 
"thinning t he chain". For an introduction to MCMC see 
iNeall (| 19931 ). for a description within a context similar to 
this one, see B07. Unfortunately, the presence of the matrix 
inverse and determinant in the likelihood function (Equa- 
tion [5| limits this algorithm to time series with less than ~ 
1500 points: even if the Cholesky decomposition is used to 
calculate det(C) and C~^y this still involves a calculation 
that takes time proportional to A*'^, where A'^ is the num- 
ber of points in the time series. For longer time series, other 
approaches are necessary; alternatively, approximations to 
Equation [5] may be possible, but are beyond the scope of 
this paper. Ad ditional effi ciency can be obtained by using 
slice sampling (|Nealll2003l ) rather than Metropolis for the 
moving of frequencies and amplitudes. 



5 SIMULATED DATA 

In this section, we demonstrate the use of our model on 
simulated data. To illustrate the method we show its out- 
put alongside output from other methods. We start with a 
simple case of a time series containing one mode and subse- 
quently expand to several modes. We generated a long time 
series by numerically solving the ordinary differential equa- 
tion [3] for a single mode of frequency 100 /xHz and damp- 
ing timescale (mode lifetime) 10^ s — 1.1574 days or 10 
oscillation periods. The amplitude of the signal was scaled 
to a standard deviation of 2 ms~^ and then evaluated at 
433 points in time, simulating 8 hours of nightly observa- 
tions over an observing period of a month. In fact, the time 
stamps w ere the same as those from the ^ Hydrae data ob- 
served bv lFrandsen et al] (|2002l ll. Thus, the simulated data 
has the same window function as the actual ^ Hydrae data. 
Measurement error was simulated by adding noise from a 
Gaussian distribution with a standard deviation of 2.5 m 
s-i. 

The results obtained from analysing this data are dis- 
played in Figure [3l The top panel shows the standard peri- 
odogram. In the 2nd panel, the results from Bayesian sine- 
wave fitting (B07) are shown. The B07 method is closely 
related to the iterative sine-wave fitting algorithm CLEAN, 
except that all frequencies are fitted simultaneously, and 
quantitative uncertainties are easily obtained. The MCMC 
approach of the B07 method actually returns a sample of 
fitted models, not a single one, however the results can be 
conveniently summarised by accumulating all detected fre- 
quencies into a single container, and then plotting a his- 
togram of the frequencies. This histogram is what appears 
in the 2nd panel of Figure [3] In the 3rd panel, a similar his- 
togram is plotted of the output from running MCMC with 
the Gaussian Process likelihood introduced in this paper. 
Finally, the "amplitude-weighted" lower panel is a similar 
histogram, however in this case, before binning, each de- 
tected frequency is given a weight proportional to its am- 
plitude. Thus, the 3rd panel indicates our confidence in the 
existence of a peak, while the 4th panel illustrates the esti- 
mated amplitude of each peak, and can be considered our 
version of an amplitude spectrum. 

The MCMC run with the Gaussian Process likelihood 
took about one hour to complete (although an MCMC run 
never really finishes, just becomes more and more useful the 
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longer it runs) on a modern PC with a 2 GHz dual core pro- 
cessor, compared to 15 minutes for the Bayesian sine-wave 
fitting and seconds for computing the periodogram. Note 
that the periodogram wou ld need further processing, such as 
Loren tzian profile fitting l|Gruberbauer. Kallinger. fc WeissI 
I2OO8I ). to obtain a posterior distribution for the frequency 
of the mode. Doing this would result in a posterior simi- 
lar to the 3rd panel of Figure [S] but with a slightly larger 
uncertainty due to the fact that the periodogram is not a 
sufficient statistic. 

The approach outlined in this paper clearly identifies 
the presence of a mode with frequency 100.38 ± 0.48 /iHz. 
This is possible because the Gaussian process model takes 
into account at the outset the fact that the predicted signal 
due to a mode is not a pure sinusoid. Lacking this informa- 
tion, the sine-wave fitting approach is forced to introduce 
many peaks in order to explain the data (Figure |3]). Note 
that the alias peaks at 90 and 110 ^Hz are automatically 
removed by the Gaussian process model. They are only par- 
tially removed by the sinewave fitting, but would have been 
completely removed had the signal been truly sinusoidal. 

A further test of this method was done by testing it 
on simulated data containing many modes. Specifically, we 
generated simulated data from a star with 11 modes with 
frequencies ranging from 50-150 /xHz in steps of 10 /xHz. 
The time series contained 433 data points at the ^ Hydrae 
times, as above. The results from analysing this simulated 
data set are shown in Figure 21 While this result is less im- 
pressive than the single mode case, the algorithm has still 
successfully identified most of the input frequencies. There 
are some anomalies, such as the merging of the peaks at 
50 and 60 /iHz, and the upward shift of the 80 /^Hz mode. 
Note that the uncertainty about each frequency can be read 
off the width of the peaks in the bottom panel of Figure H) 
Whilst our method provides cleaner results than the raw 
periodogram, it is clearly not perfect. When interpreting re- 
sults from this method, the summary plots like those shown 
in Figure [4] may be used as a guide, but the full output of 
the MGMC sample should be considered when the results 
are critical. 

The mode lifetime r can be measured using our anal- 
ysis, as it is just another parameter that gets estimated by 
the MCMC. The posterior distribution for the mode life- 
time is simply a histogram of the r values encountered by 
the MGMC chain. For the multiple-mode simulated data, 
this distribution is shown in Figure O The true input value 
of 10^ s = 1.1574 days is recovered, albeit with a large un- 
certainty, which is unsurprising given the time series is only 
433 points in size. The distribution is asymmetric, largely 
due to our choice of a 1/r prior. Thus, conventional error 
bars are inappropriate. An alternative statement of uncer- 
tainty is the symmetric 95% credible interval for the mode 
lifetime, which is [0.58, 2.10] days. 



6 THE LARGE SEPARATION AND MODE 
LIFETIME OF HYDRAE 

We now turn to the real observations of ^ H ydrae. The first 
preliminary analysis of the data presented bv lFrandsen et al] 
(|2002l ) showed strong evidence for solar-like oscillations 
based on the amplitude, the frequency range, and the fre- 
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Figure 3. Results from analysing the single mode simulated data 
with various methods. Bayesian sine-wave fitting (2nd panel) can 
only explain the time series data by introducing many peaks, 
whereas the actual simulation contained only a single input fre- 
quency. Using the more realistic Gaussian Process likelihood, we 
find that a single frequency can explain the data (lower two pan- 
els). 



quency separation of the extracted modes, which all agreed 
with theoretical predictions. They assumed that the mode 
lifetime was re latively long, in accord with the theoretical 
calculations bv lHoudek fc Gouehl (|2002l ') (r ~ 17 days), and 
hence their analysis relied on the conventional power spec- 
trum and iterative sine- wave fitting (GLEAN). The value 
they found for the dominant frequency separation was 7.1 
/iHz found between the strongest modes and 6.8 /iHz found 
from an autocorrelation of the power spectrum in the re- 
gion of excess power. Subsequent studies of th e same data 
including extensive simulations bv lStello et al] (|2004l ) indi- 
cated that the mode lifetime was significantly shorter than 
the theoretical value (r 2 days). This result was further 
confirmed bv iStello et al.l ([200 61 using an independent ap- 
pr oach which also confirm ed the frequency separation found 
bv lFrandsen et al] (|2002l ). but they showed that the preci- 
sion by which the frequency separation could be established 
from the data was low due to the short mode lifetime. In 
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Figure 4. Results from analysing simulated data witii input fre- 
quencies from 50-150 fiRz in steps of 10 fiHz. Most of the modes 
are successfully identified, although some are spuriously shifted. 
The dashed lines indicate the true input frequencies. 
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Figure 5. The posterior distribution for mode lifetime t, given 
the simulated data set with 11 modes. The uncertainty is quite 
large, but comfortably contains the true input value 10^ s = 
1.1574 days. 



panel of Figure [7] - the area under the curve over any fre- 
quency range is proportional to the probabihty that a mode 
exists within that range, yet even the peaks in this plot are 
only a factor of ~ 2-3 higher than the background. There 
is some suggestion of a regular pattern to the peaks. To 
measure the large frequency separation, we took the power 
spectrum of the lower panel of Figure [7] (a full Bayesian es- 
timate of the large frequency separation, as done in B07, is 
prohibitive in this case) in order to search for periodicities. 
This power spectrum is displayed in Figure |8] and shows at 
least three possible periodicities in the frequency pattern: 
one at 6.3 /iHz, another at 9.6 /iHz and a third peak at 
19.2 jjJiz (although this is simply a doubling of the 9.6 /iHz 
peak). Usually, if / = and I = 1 modes are detected, the 
dominant separation of modes is half of the large frequency 
separation. This would imply that the large separation of 
^ Hydrae is 12.6, 19.2 or 38.4 /iHz. 

However, from the classical stellar parameters (luminos- 
ity, mass, and effective temperature), the estimated large 
separation is 7.0 /xHz wit h about 10% uncertainty, us ing the 
solar scaling presented in lKieldsen fc Beddind l| 19951 ). Also, 
a stellar pulsation model that goes through the star's po- 
sition in th e H-R diagram gives an average large spacing 
of 7.2 /iHz (|Frandsen et al.l 120021 '). Hence, the most plausi- 
ble solution consistent with stellar astrophysics is that 6.3 
/iHz is the large separation, not half of the large separa- 
tion. Thus, we conclude that radial modes contributed most 
of the signal, non-radial modes may b e excited with am - 
plitudes below the detection threshold jStello et al.ll2006l ). 
Although we have not formally modelled the uncertainty in 
the large separation, inspection of Figure [6] shows that the 
uncertainty must be large. This uncertainty may be reduced 
with further observations, or perhaps by taking theoretical 
models of the star into account as prior information (B07). 
Our analysis used a uniform prior distribution for the fre- 
quencies, but if a large ensemble of plausible stellar models is 
produced, this would significantly reduce the range of possi- 
ble frequency patterns and significantly improve the quality 
of the data analysis. Performing such an analysis is beyond 
the scope of this paper. 

The posterior distribution for the mode lifetime of 
^ Hydrae is shown in Figured The mode lifetime is found to 
be very short, of order 1 day, albeit with a fairly large uncer- 
tainty. An estimate with 1-a error bars is logj^Q(r/l day) = 
0.03±0.21, and we find that r lies between 0.41 and 2.65 days 
with 95 % posterior p r obabi lity. This compares well with the 
result of IStello et al.l (|2006| ) who estimated r to be 2 days, 
but also with a large uncertainty. This mode lifetime remains 
much shorter than the theoretical predictions of 15-20 days 
(Houdek fc Gough,2002. ). 



this section we will apply our new Gaussian Process method 
to the ^Hydrae observations and compare our results with 
those found by the previous studies. 

The results from running our code on the ^ Hydrae data 
are displayed in Figures [6] and [T] The diversity of the mod- 
els in Figure |6] indicates that the uncertainties are quite 
large, and only a few modes are securely detect ed; this re- 
sult agrees with the analysis of lStello et al.l l|2006l ) . The large 
uncertainty about the frequencies is confirmed by the lower 



7 CONCLUSIONS 

In this paper, we have described a new Bayesian method for 
inferring the frequencies and amplitudes (with uncertain- 
ties on everything, including the number of modes present) 
of stellar oscillation modes from time series observations of 
the radial velocity or the intensity of the star. The method 
includes a Gaussian Process likelihood, which allows us to 
take into account the fact that the predicted signature of 
an oscillation mode is not exactly sinusoidal. Exactly how 
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Figure 6. A sample of models from the posterior distribution for 
the oscillation frequencies and amplitudes of ^ Hydrae. Clearly, 
the data and prior information do not uniquely determine the cor- 
rect model. However, any question about the frequencies present 
can be answered probabilistically by calculating the fraction of 
the output models that have the property that is being tested 
for. The full sample is much larger than the nine models shown 
here. 



Conventional Power Spectrum 




20 40 60 80 100 120 140 160 180 200 

Amplitude-Weighted 




20 40 60 80 100 120 140 160 180 

Probabiiity-Weighted 




20 40 60 80 100 120 140 160 18 

Frequency (p. Hz) 



Figure 7. Summarised results for the frequency spectrum of 
^ Hydrae. 
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Figure 8. Power spectrum of the estimated frequency spectrum 
(lower panel of Figure (T]! of § Hydrae. A regular pattern to the 
frequencies of the modes should show up as a peak in this plot. 
There are high peaks at 6.3 ^Hz, 9.6 /iHz and 19.1 /iHz. 
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Figure 9. The posterior distribution for mode lifetime t for 
5 Hydrae. The probability distribution for t is asymmetric, but 
for log(r) it is approximately Gaussian. Hence, we estimate 
logj^Q(r/l day) = 0.03 it 0.21 (1-sigma error bars), r lies between 
0.41 and 2.65 days with 95% posterior probability. 



non-sinusoidal the oscillation signals are, is described by the 
mode lifetime, which is also estimated from the data, along 
with a measurement of the uncertainty in this value. 

The method was implemented using a Markov Chain 
Monte Carlo algorithm and applied to two simulated data 
sets. As expected, the method removed the extra peaks 
caused by aliasing and the finite mode lifetime. We spec- 
ulate that this method is, at the very least, comparable to 
the results obtained by fitting the power spectrum, but is 
more straightforward to interpret. Our method also avoids 
any concerns about information loss due to the fact that the 
power spectrum is not a sufficient statistic; whether this is 
of significant practical importance depends on the sampling 
of the time series, and the mode lifetime. For well-sampled 
time series and long mode lifetimes, information loss is not 
an issue. 

Unfortunately, the method presented in this paper is 
computationally intensive due to the presence of a matrix 
inverse and determinant in the likelihood function. This lim- 
its its practical use to small time series with less than about 
1500 points. For time series with 1500-15000 points, the 
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Bayesian sine wave fitting approach (B07) is recommended, 
and for tliose witli more tlian 15000 points, neitlier is com- 
putationally feasible; periodogram-based analysis is clearly 
the best choice here. 

Applying the method to radial velocity data of the red 
giant ^ Hydrae, we found that the mode lifetime lies be- 
tween 0.41 and 2.65 days with 95% posterior probability. 
The large frequency separation was estimated to be either 
6.3 fiHz or 9.6 /iHz, with the former being the most plausible 
given the star's position in the Hertzsprung-Russell diagram. 
C-|--|- programs implementing the methods described in this 
paper (both the Gaussian Process and the sinewave fitting 
versions) are available upon request from B. J. Brewer. 
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